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Abstracts  We-  propose  a  new  method  for  the  solution  of  the  wide  angle  wave  equation  in 
three  dimensions.  In  contrast  wjtn  standard  techniques,  our  approach  requires  only  solu¬ 
tions  of  successive  tridiagonal  systems  in  the  resulting  finite  difference  parabolic  equations. 
The  method  is  basgd-oft  a  simple  approximation  to  the  square-root  operator  written  for¬ 
mally  as  tfT+X  +  Y  where  X  is  a  partial  differential  operator  with  respect  to  the  depth 
z  and  Y  is  a  partial  differential  operator  with  respect  to  the  azimuthal  angle  Q.  W^exploits 
the  fact  that  the  partial  derivative  term  Y  with  respect  to  the  azimuthal  angle  is  small, 
but  not  negligible,  as  compared  with  other  terms.  It  is  then  natural  to  replace  the  square- 
root  operator  by  an  expansion  which  is  of  order  2  with  respect  to  the  X  operator,  and 
of  order  1  with  respect  to  the  Y  operator.  An  important  feature  of  this  approach  is  that 
it  is  then  possible  to  derive  a  rational  function  approximation  to  the  exponential  of  the 
square-root  operator  which  has  the  property  of  being  stable,  and  accurate.  Moreover,  the 
approximation  decouples  naturally  as  a  product  of  a  (1, 1)  rational  function  of  „Y  times  a 
(1, 1)  rational  function  of  Y.  As  a  consequence,  this  will  result  in  a  solution  technique  that 
requires  only  two  tridiagonal  system  solutions  per  step,  namely  one  for  the  X  operator 
and  one  for  the  Y  operator.  Numerical  examples  are  reported  that  show  the  wide  angle 
capability  of  this  method.  ( 


O 

* 


u, 

dc 


An  efficient  method  for  solving  the 
three-dimensional  wide  angle  wave  equation 


Ding  Lee,  Youcef  Saad  and  Martin  H.  Schultz 

Research  Report  YALEU/DCS/RR-463 
October  1986 


?h 


has  been  approved-)^ 

di  lr.laase  and  sale;  its 

.  J;;  '  n  n  unlimited. 


DT1G 

tf&eLECTE| 

JUN  0  91987  J  g 

fp 

E 


This  work  was  supported  in  part  by  the  Office  of  Naval  Research  under  contracts  N00014- 
82-K-0184,  N00014-85-WR-240G8,  N00014-85-WR-24268,  N00014-86-WR-24233,  ancl  in 
part  by  Naval  Underwater  Systems  Center  Independent  Research  project  A65020. 


1.  Introduction 


The  wide  angle  three-dimensional  parabolic  approximation  technique  developed  by 
Siegmann,  Kriegsmann,  and  Lee  [5]  has  been  proven  capable  of  handling  wide  angle  prop¬ 
agation  in  the  vertical  plane.  This  technique  is  based  on  a  pseudo-differential  three-  di¬ 
mensional  parabolic  wave  equation  of  which  the  3-D  parabolic  approximation  introduced 
by  Tappert  [7]  is  a  special  case.  Methods  for  solving  this  equation  have  been  developed 
by  Baer  and  Perkins  [l],  using  the  fast  Fourier  transform,  but  are  only  applicable  to  small 
angles  of  propagation.  Moreover,  these  methods  do  not  extend  easily  to  equations  with 
variable  coefficients  and  do  not  handle  rigid  boundary  conditions. 

In  order  to  be  able  to  handle  wide  angle  propagation  in  the  general  variable  coeffi¬ 
cient  case  and  to  easily  treat  rigid  boundary  conditions,  we  choose  the  three-dimensional 
parabolic  wave  equation  as  the  representative  equation.  The  main  contribution  of  this  pa¬ 
per  is  a  method  for  efficiently  solving  this  equation  for  wide  angle  propagations.  A  solution 
technique  for  the  wide  angle  3-D  equation  has  been  previously  developed  by  Schultz,  Lee, 
and  Jackson  [4]  using  the  Crank-Nicolson  scheme  in  conjunction  with  a  preconditioned 
conjugate  gradient  method.  Due  to  the  resulting  properties  of  the  discretized  finite  differ¬ 
ence  equations,  the  operator  is  neither  Hermitian  nor  positive  real.  Therefore  no  effective 
preconditioners  are  known  for  this  case  and  the  solution  adopted  by  Schultz,  Lee,  and 
Jackson  [4]  was  to  precondition  the  normal  equations.  This  squares  the  condition  number 
of  the  initial  matrix  but  gives  satisfactory  results  as  far  as  accuracy  is  concerned. 

In  this  paper  we  propose  a  new  method  to  solve  the  same  3-D  wide  angle  parabolic 
approximation.  What  makes  our  technique  so  attractive  as  compared  to  other  techniques 
is  that  each  integration  step  requires  solving  only  two  successive  tridiagonal  systems.  The 
method  resembles  alternating  direction  schemes  but  its  foundation  and  analysis  are  differ¬ 


ent.  The  main  goal  of  this  paper  is  to  derive  this  method  and  to  discuss  its  validity  and 
accuracy.  The  theory  is  then  verified  by  performing  two  numerical  tests,  an  azimuthally  r 
independent  case  and  azimuthally  dependent  one.  The  first  example  is  for  testing  the 
accuracy  of  the  method  and  for  determining  how  wide  an  angle  it  can  accomodate.  The 
second  example  is  for  verifying  whether  angular  dependencies  are  well  handled  and  for  “ 
comparing  the  speed  of  our  method  with  the  speed  of  other  methods.  This  last  test  shows  ^ 
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an  example  where  the  new  method  is  orders  of  magnitude  faster  than  existing  competing 
techniques. 


2.  Background 

The  standard  wide  angle  3-D  wave  equation  was  developed  using  the  classical  formu¬ 
lation  of  the  Helmholtz  equation  in  three  dimensions,  in  cylindrical  coordinates  ( r,0,z ): 


d2p  1  dp  1  d2p  d2p  22  ^ 

1?  +  rTr  +  +  d?  +  *°n  P~°- 


(2.1) 


In  the  abo/e  equation  p  represents  the  acoustic  pressure,  ko  =  w/co,  where  co  is  a  reference 
sound  speed,  w  =  2rrf,  in  which  /  is  the  frequency  of  the  signal  and  finally  n  =  n(r,  9 ,  z )  = 
co/c(r,  9 ,  z)  is  the  index  of  refraction,  in  which  c(r,  9 ,  z)  is  the  sound  speed. 

We  make  a  standard  transformation  of  the  above  equation  by  writing  the  pressure  in 
the  form  [7]: 

p{r,9,z)  =  u(r,9,z)v{r) 


where  the  factor  v(r )  represents  a  rapidly  varying  portion  of  the  pressure  and  u(r,  9,  z)  is 
its  modulation,  a  slowly  varying  function  with  respect  to  range.  After  neglecting  small 
terms,  making  use  of  the  far-field  approximation  [k^r  »  l),  and  rearranging  the  above 
equation  we  obtain: 


d2u  _..du  1  d2u  d2u  ,  2  .\i2 


—  +  2/fco  — — h 


+  Tpy  +  (n  —  l)^o“  —  0- 


(2.2) 


dr*  '  ~"^dr  r2d9 2  dz 2 

This  new  equation  has  been  at  the  origin  of  the  very  successful  small  angle  parabolic 
approximation  technique,  which  consists  in  simply  dropping  the  second  order  derivative 
with  respect  to  r  and  integrating  the  resulting  parabolic  equation. 

It  is  convenient  to  define  the  operators: 

1  d2 


X  = 


A’q  dz2 


+  («  -  1), 


Y  = 


i  a2 

*2r2  d92' 


after  which  the  above  equation  reads  as  follows: 


d2u  du  .  o  .  __ 

+  k0  (X  +  Y)  u  =  0. 


(2.3) 

(2.4) 

(2.5) 
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The  standard  wide  angle  PE  technique  [5]  starts  by  approximately  factoring  the  above 
operator  as  the  product 


-r-  +  iko  -  ikoQ 
or 


ill 


4-  ik0  +  ikoQ 


in  which 


Q  =  y/1  +  X  +  Y. 

Then  the  operator  Q  is  approximated  by  a  rational  function  in  the  form 

^  I  +  P1X  +  P2Y 
V  ~  1  +  qiX  +  q2Y 

which  yields  the  wide  angle  ‘parabolic’  equation 


(2.6) 


(2.7) 


(2.8) 


f  1  +  Pl^  +  P2^\ 

“r  =  r  0 


(2.9) 


To  solve  (2.9)  Schultz,  Lee  and  Jackson  (4]  applied  the  Crank-Nicolson  scheme  in  conjunc¬ 
tion  with  a  preconditioned  conjugate  gradient  method.  The  application  of  Crank-Nicolson 
reduces  (2.9)  to  a  sequence  of  systems  of  difference  equations  of  the  form 


(/  -  ^ArL)un+1  =  (/  +  ^A rL)un, 

&  mt 


where 


L  =  —iko  tfco 


1  -I-  piX  +  P2Y 


(2.10) 


(2.11) 


1  +  <71 X  -f-  <72^ 

By  multiplying  both  members  of  (2.10)  by  the  denominator  of  (2.11),  we  obtain  a  marching 
process,  in  which  a  large  block-tridiagonal  linear  system  must  be  solved  at  each  step. 


3.  The  new  approach 

Our  approach  starts  with  the  equation  (2.5)  which  is  formally  considered  as  an  ordi¬ 
nary  differential  equation  with  respect  to  the  variable  r.  For  convenience  the  variables  z 
and  6  will  therefore  be  dropped  out  in  the  remainder  of  the  paper:  u(r)  stands  for  u(r,  z,  9). 
Locally,  its  formal  solution  has  the  form 

u(r  +  Ar)  =  e-**oAre«*0Arv/I+X+F  +  e-ik0Are-ik0^T+Y+Yu-^ 
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where  u+(r)  and  u~(r)  are  some  initial  conditions  at  the  range  r.  The  first  term  in  the 
above  solution  is  the  outgoing  wave  and  the  second  is  the  incoming  wave.  In  this  paper 
we  will  neglect  back-scattering  and  therefore  the  second  term  will  be  dropped  to  yield  the 
local  solution 

u(r  +  A  r)  =  e~ses^l+x+Yu(r)  (3.1) 


in  which  we  have  set 


6  =  t'fcoAr. 


Note  that  this  is  also  a  local  solution  of  the  one  way  wave  equation 


ur  =  (-iko  -(-  ifco'/l  +  X  +  Y)u 


which  is  obtained  by  neglecting  the  second  factor  in  (2.6). 

The  approach  taken  in  this  paper  consists  of  approximating  the  term  eWi+X+Y  jn 
a  convenient  and  accurate  manner.  An  easy  way  in  which  this  can  be  done  is  to  use  the 
approximation 

Vl  +  X  +  Y  »  1  +  hx  +  Y) 

& 

which  yields  the  standard  three-dimensional  narrow  angle  parabolic  equation: 


ur 


i  <92 

2k^d^  + 


i  d 2  \ 

2J t0r2  ae 2 ) 


u  =  Lu. 


This  equation  was  solved  by  Baer  and  Perkins  [1,3]  using  a  split-step  Fourier  algorithm. 
However  this  equation  accurately  represents  only  narrow  angle  propagation. 

To  accomodate  wide  angle  propagation,  we  consider  the  higher  order  approximation 


Vl  +  X  +  Y  *s  1  +  -X  -  - X 2  +  -Y. 

2  8  2 


i  +  \x  -  i y2  +  irh  u , 
2  8  2  \) 


The  corresponding  approximation  to  the  wave  equation  becomes 

«r  =  ^—ik o  -f  iko 

and  formula  (3.1)  becomes: 

u(r  +  Ar)  =  e-SeS^hX-hX2+hY)uir) 


(3.2) 


(3.3) 
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(3.4) 


Assuming  that  n(r,  6 ,  z)  varies  slowly  with  respect  to  0 ,  the  operators  X  and  Y  are 
nearly  commutative,  and  equation  (3.3)  yields 

u{r  +  Ar)  =  e"V(1+SX“8*2)e5ru(r)  (3.5) 


In  the  following  we  seek  an  approximation  of  the  term 


(3.6) 


Here,  the  function  G  and  its  approximations  should  be  regarded  as  functions  of  the  real 
variable  Ar  while  6  is  an  independent  parameter. 

The  Taylor  series  expansion  of  (3.6)  about  X  =  0  is 

1  /  <52  6' 


G(6 ,  X)  =  e6  [l  +  6-X  +  i  (i-  -  0  X2]  +  0(X3). 


(3.7) 


We  seek  an  approximation  to  the  function  G  in  the  form 

G^”‘llTTW' 


(3.8) 


where  p  is  a  complex  number  to  be  determined. 

Writing  that  the  expansion  (3.7)  is  equal  to  the  right  hand  side  of  (3.8)  up  to  0(X3) 
we  get  the  equation 

1  +  *V  +  i^-^  Y2  =  i±M 

+  2‘  +2!\4  4 /  ^  1  4-  pX 

from  which  it  is  easy  to  obtain  p: 


I  6 
P=4  +  4 


A  similar  development  for  the  term 


H{6,Y)  =  e? 


-  JjY 


leads  to  the  approximation 


with 


H(6,Y) 


1  +  qY 
1  +  qY 


(3.9) 


q=4- 
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(3.10) 


Therefore  we  get  the  final  expression 

i + a  ±  $)* 


u(r  +  A  r)  = 
This  can  be  rewritten  as 


+  (?  “  I 

u(r  +  A r)  =  L  u(r) 


fl+fl'l 

^4 

'-ohf 

i 

H 

_ i 

u(r). 


where 

L  =  LxLy*LY,  (3*11) 

in  which 

Lx  =  1  +  (j  +  j)  X,  LY  =  1  +  jV,  (3.12) 

and  where  X-*  stands  for  the  inverse  of  X*,  the  adjoint  of  X.  Equation  (3.10)  can  formally 
regarded  as  an  explicit  marching  scheme.  It  can  easily  be  seen  that  the  two  operators  in 
the  denominator  of  (3. 10)  are  nonsingular  because  <5  is  purely  imaginary  and  X  and  Y  are 
both  self-adjoint.  In  the  next  section  we  analyze  the  accuracy  and  stability  of  this  scheme 
and  in  Section  5,  we  will  see  how  to  discretize  it  and  use  it  numerically. 


4.  Theoretical  aspects 


4.1.  Stability 

The  operators  X  and  Y  defined  earlier  are  self-adjoint  and  therefore  their  correspond¬ 
ing  eigenvalues  are  real.  Since  the  numerator  and  denominator  of  each  term  between 
brackets  in  equation  (3.10)  are  the  conjugate  of  each  other,  a  modal  expansion  of  (3.10) 
shows  that  with  respect  to  each  mode,  the  error  induced  by  the  marching  scheme  will  not 
increase  exponentially,  i.e.,  the  scheme  (3.10)  is  stable.  After  discretization,  the  eigenval¬ 
ues  of  X  and  Y  will  remain  real  provided  boundary  conditions  are  properly  handled  and 
the  discretizations  in  the  numerators  and  denominators  are  calculated  at  the  same  range  r 
(see  next  section).  Under  these  conditions  the  scheme  is  stable.  Note  that  the  second  part 
of  (3.10)  is  the  usual  Crank- Nicolson  approximation  applied  here  to  the  term  e^/2)y, 

4.2.  Accuracy 

To  analyze  the  local  error  of  the  integration  scheme  (3.10)  we  must  attempt  to  find 
an  estimate  of  the  difference  between  the  operator 

e-6e6y/l+X+V  (4  .i) 
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and  the  operator  in  the  right  hand  side  of  (3.10).  In  the  following  we  consider  that  X  and 
Y  are  two  independent  real  variables.  On  the  one  hand,  we  find  after  some  calculation 
that  the  second  order  Taylor  expansion  of  the  operator  (4.1)  is 


e-6e6y/l+X+Y  =  1  +  +  Y)  +  1(6  -  1) 

2  o 


X2  +  2  XY  +  K2]  +  0(||(X,  F)||3)  (4.2) 


On  the  other  hand,  the  second  order  Taylor  expansion  of  the  operator  of  the  right  hand 
side  of  (3.10)  is  given  by 


'i  +  (J  +  £)*' 

'l  +  f  Y 

r 

,i  +  (i  -  i)x. 

L 

1 + !* + -  w2 + -■ 


6  62  o 

1  +  2K  +  Ty  +  -• 


=  1  +  6-(X  +  Y)  +  J XY  +  -  1)X2  +  J Y*  +  .. 

The  larger  terms  in  the  difference  between  the  two  (4.2)  and  (4.3)  are 


(4.3) 


,—6  6\/l+X+Y  _ 


'l  +  {y' 

L*-H 

=  -{xy  -  jF2  +  o(|K.r,  y)H3). 


Thus,  one  can  expect  good  accuracy  when  Y  is  very  small  and  X  is  small.  The  above 
error  is  better  than  an  error  of  the  form  0(X2),  because  the  error  expression  is  the  product 
of  three  small  terms  namely  6,Y,  and  X.  Moreover,  we  have  assumed  that  the  term  Y  is 
much  smaller  than  the  term  A'. 

When  using  a  marching  scheme,  an  upper  bound  for  the  global  error  at  some  point 
can  be  derived  from  the  above  local  truncation  error  by  expressing  the  error 


ej+1  =  -  J+K 

where  uJh+1  represents  the  computed  solution  at  step  j  ai_d  u3+l  the  exact  solution  at  the 
same  point.  On  the  one  hand,  we  have 

uh+l  =  Lhui 


where  L /,  is  the  discretization  of  the  operator  L  as  defined  by  (3.11).  On  the  other  hand 

u;+1  =  LhU}  +  tj 
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where  Cj  is  the  truncation  error  incurred  at  step  j.  Hence, 


ei+1  =  Li>ei  +  e) 

A 

and  since  the  operator  L/,  is  unitary,  we  have 


i 


K+,ii<EiM- 

i=0 


Since  the  number  of  steps  in  range  is  0(5  1),  the  above  analysis  shows  that  the  norm  of 
the  global  error  ||eji+1||  is  0(A'2). 


5.  Finite  Difference  Solution 

To  employ  formula  (3.10)  numerically,  we  must  discretize  the  operators  X  and  Y  by 
central  differences  and  replace  the  corresponding  operators  by  their  discrete  analogues. 
However,  we  first  put  equation  (3.10)  in  the  form 


uJ+1  = 


,1  6, 

6  1 

-l 

r  5 

l  +  (-  +  -)X 

[i+irJ 

(5.1) 


There  are  several  way  of  rewriting  (5.1)  exploiting  the  commutativity  of  the  operators 
L\  and  LJr*  and  of  Ly  and  Ly*.  The  near  commutativity  of  X  and  Y  can  also  be  exploited 
to  derive  alternative  formulae.  For  example  we  may  consider  the  scheme 


1  +  (j  -  j)X 

-  V 

u>+1  = 

1  +  ( j  +  j)X 

1  +  jY 

4  4 

4 

4  4 

4 

Although  not  obvious  at  first,  the  above  scheme  is  also  unconditionally  stable.  The  reason 
for  this  is  that  the  corresponding  operator  L  =  Ly*  LJ?  Lx  Ly  is  also  unitary  as  is  readily 
seen  by  forming  LL*  which  is  found  to  be  the  identity  operator.  The  question  as  to  which 
of  the  various  schemes  is  tc  be  preferred  is  certainly  worth  further  investigation  but  we 
will  not  pursue  it  in  the  present  paper. 

Let  us  denote  by  A  the  finite  difference  approximation  of  the  operator  Lx  =  I  +  — 

|)AT  and  by  B  that  of  Ly  =  I  +  \Y .  Both  matrices  are  tridiagonal  with  the  structure 
indicated  in  Table  1. 
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Super  diagonal 

Diagonal 

Sub  diagonal 

A 

1  +  (?  “  ?)("2  “  x)  “  2(?  ” 

(i-i)j 

B 

6  1  1 

t  ,  <5  l  l 

1  + 

6  1  1 

3  k y  (A 0)s 

Table  1:  Coefficients  of  the  two  tridiagonal  matrices  A 
and  B 


When  solving  tridiagonal  systems  with  the  matrices  A  and  B  it  is  of  interest  to  know 
whether  these  matrices  are  diagonally  dominant  or  not.  While  the  matrix  B  is  always 
diagonally  dominant,  the  situation  is  more  complicated  for  A.  In  the  simple  case  where 
n(r,6,z)  =  1,  the  matrix  is  conditionally,  i.e.,  for  h  >  1  /ka,  diagonally  dominant  in  the 
sense  that  the  modulus  of  the  diagonal  term  is  greater  than  or  equal  to  the  sum  of  the 
modulii  of  the  off-diagonal  terms  in  the  same  row.  The  more  general  case  where  n  is 
arbitrary  is  not  easy  to  analyze. 

The  scheme  (5.1)  becomes 

u>+1  =  A—AB—Bui.  (5.3) 

Note  that  we  evaluate  the  matrices  A  and  B  at  mid  distance  between  and  uJ ,  i.e.,  at 
range  r  +  Ar/2.  This  is  in  order  to  ensure  that  the  operators  A*  and  A,  as  well  as  B *  and 
B ,  form  two  pairs  of  operators  that  are  conjugate  of  each  other.  This  choice  will  guarantee 
stability  as  was  seen  in  Section  4.1 

To  perform  one  step  of  (5.3)  we  must  start  by  computing  iuJ  =  Bu}  and  solve  the 
tridiagonal  system 

£  V+1  =  w1'. 

Then  we  compute  w 1  =  A~1A*uJ+l  and  solve  the  tridiagonal  system 

A*u3+l  =  AwJ.  (5.4) 

Thus  there  are  two  multiplications  of  a  tridigonal  system  by  a  vector  and  two  tridiagonal 
systems  to  solve  at  every  step. 
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6.  Test  examples 

The  numerical  scheme  (5.3)  has  been  implemented  into  a  research  computer  code 
which  is  used  to  predict  wave  fields  at  required  ranges.  We  used  the  marching  scheme  (5.2) 
instead  of  the  original  scheme  (5.1).  The  resulting  values  are  compared  against  a  known 
solution  to  check  the  validity  as  well  as  the  accuracy  of  our  scheme.  All  computations 
are  made  on  the  VAX-11/780  computer  using  single  precision  complex  arithmetic.  We 
present  two  examples.  The  first  one  is  an  azimuthally  independent  case  and  the  second 
is  an  azimuthally  dependent  one.  The  input  parameters  for  both  examples  are  shown  in 
Table  2  for  convenience. 

6.1.  An  azimuthally  independent  case 

To  start  the  computation,  the  initial  field  is  taken  from  the  following  formula  borrowed 
from  [6]: 


where  ay  satisfies 


We  are  concerned  with  the  propagating  mod°s  namely  those  for  which  ay  remains  real. 

In  the  computation  the  sector  is  divided  into  10  portions.  The  two  extreme  sectors 
represent  sector  boundaries  where  the  solution  is  supplied  from  the  exact  solution.  The 
inner  sectors,  because  of  the  azimuthal  independence  have  the  same  initial  values.  It  is 
expected  that  at  different  sectors  at  the  same  depth,  the  computed  wave  fields  should  be 
identical.  We  examined  this  hypothesis.  Another  verification  we  did  was  to  look  at  the 
size  of  the  angle  of  propagation.  To  simulate  this  we  took  the  mode  index  j  to  be  9  so  that 
we  could  obtain  an  angle  of  propagation  of  around  52°.  Our  method  could  handle  such  a 
wide  angle  without  any  major  difficulty. 

Our  results  are  summarized  in  the  two  tables  3  and  4.  The  first  table  shows  the 
accuracy  achieved  and  the  second  shows  the  angle  of  propagation  for  the  case  of  8  modes 
and  9  modes.  In  the  tables  the  values  appearing  in  the  first  row  represent  the  calculated 
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Input  parameters 

Problem  1 

Problem  2 

Source 

300  m 

10  m 

Initial  range 

100  m 

10  m 

Source  frequency 

20  Hz 

20  Hz 

Bottom  depth 

400  m 

20  m 

Sound  speed 

1500  m/s 

1500  m 

Reference  sound  speed 

1500  m/s 

1500  m/s 

Receiver  depths 

156  m,  312  m 

5  m 

Propagating  sector 

— 2^,  +2ij 

-5,  5 

Depth  increment 

4  m 

0.2  m 

Range  step  size 

1  m 

0.001  m,  0.25  m 

Angular  increment 

0.5° 

1° 

Maximum  range 

1  km 

10.5  m 

Surface  condition 

pressure  release 

Dirichlet 

Bottom  condition 

Rigid 

Dirichlet 

Size  of  matrices  A  and  B 

1000 

1000 

Table  2:  Parameters  for  the  two  test  problems. 


e 

zr{m) 

—  2  degree 

+  2  degree 

156 

-0.59261E-04  0.20995E-04 
-0.59224E-04  0.20954E-04 

-0.59259E-04  0.20995E-04 
-0.59224E-04  0.20954E-04 

312 

-0.96924E-04  0.34283E-04 
-0.96908E-04  0.34287E-04 

-0.96924E-04  0.34285E-04 
-0.96908E-04  0.34287E-04 

Table  3:  Wave  field  results  at  1  km  range:  accuracy 
test. 


values  by  our  new  method,  the  values  in  the  second  row  are  the  exact  solutions.  The  first 
columns  are  the  real  parts  and  the  second  columns  are  the  imaginary  parts. 

6.2.  An  azimuthally  dependent  case 

This  second  example  deals  with  a  »ow  frequency  propagation  in  shallow  water.  To 
construct  an  azimuthally  dependent  case,  we  modified  a  reference  solution  tested  by  Chan, 
Shen  and  Lee  [2]  and  used  the  same  exact  input  parameters  to  derive  a  system  of  equa¬ 
tions  with  the  same  dimension  in  order  to  compare  the  computation  speed.  An  exact 
solution,  after  the  modification,  to  the  wide  angle  3-dimensional  wave  equation  (3.3)  can 


11 


Mode  j 

Angle  size 
(Degrees) 

Results 

Complex  values 

Results 

dB 

6 

31.03° 

(0.1374SE-04, -0.23722  E-05) 
(0.15624E-04, -0.12743  E-05) 

97.107 

96.095 

7 

37.54° 

(-0.96649E-05,  0.25102  E-05) 
(-0.68974E-05,  0.23300  E-05) 

100.013 

102.757 

Table  4:  Wave  field  results  at  1  km  range:  angle  of 
propagation  measurements. 


be  expressed  by 

u(r,  9 ,  z)  =  C-nVmV‘m2/(2*or)  (6i3) 

The  expression  is  used  to  generate  the  initial  field.  The  input  parameters  used  are  shown 
in  the  second  column  of  Table  2.  The  surface  condition  is  taken  to  be 

ti(r,^,0)  =  e,mVm2/(2M  (6.4) 

and  the  bottom  condition  was 

«M,*max)  =  e"n*™*eimVm2/(2*°r)  (6.5) 

The  scalar  fi  in  Equation  (6.5)  is  chosen  to  be  2k$.  The  angular  modal  number  m  is  taken 
to  be  3.  To  obtain  an  accuracy  of  10-2,  as  in  [2],  their  methods  need  to  take  a  range  step 
size  of  0.001  m. 

We  tested  two  different  step  sizes:  0.001  m  and  0.25  m.  The  experiment  with  the  first 
range  step-size  is  only  done  for  a  comparison  with  reference  [2j.  We  should  point  out  that 
such  a  small  step  size  for  the  5-point  method  of  Chan,  Shen  and  Lee  is  necessarv  because 
the  scheme  is  explicit.  In  this  computation,  we  found  that  our  method  was  approximately 
1.6  times  faster  than  the  5-point  explicit  method  of  [2],  and  17  times  faster  than  Crank- 
Nicholson  of  [4].  Note  that  the  Crank-Nicholson  of  [4]  uses  a  stable  version  of  the  conjugate 
gradient  method  applied  to  the  normal  equations,  called  Craig’s  method.  Using  the  second 
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Method 

Ar 

Relative  error 

CPU  time  (h-m-s) 

Crank  Nicolson 

0.001 

(0. 18E-01  ,-0. 12E-01) 

03-47-10 

5-Point  Explicit 

0.001 

(0.10E-01,-0.11E-01) 

00-21-35 

New  method 

0.001 

(0.26E-02,-0.1lE-02) 

00-13-12 

New  method 

0.25 

(0.22E-02,-0.1lE-01) 

00-00-09 

Table  5:  Results  for  Problem  2. 

step  size  of  0.25  m,  we  found  that  the  same  accuracy  could  be  achieved  by  our  method  as 
with  A r  =  0.001  but  the  excution  was  much  faster.  Here,  our  method  is  approximately  160 
times  faster  than  the  5-point  method  and  1600  times  faster  than  Crank-Nicolson.  Results 
are  displayed  in  Table  5. 

7.  Conclusion 

Obtaining  solutions  to  ocean  acoustic  propagation  in  three  dimensions  can  be  very 
complicated  and  computationally  expensive.  Moreover,  it  is  now  becoming  the  general 
consensus  that  two-dimensional  models  are  no  longer  sufficiently  representative.  Efficient 
methods  and  clever  implementations  for  dealing  with  three-dimensonal  wave  propagation 
are  therefore  very  important. 

The  new  method  proposed  in  this  paper  is  not  only  a  fast  and  accurate  method,  but 
also  has  the  property  of  being  as  representative  a  model  as  other  well  known  existing  3-D 
models.  Our  numerical  results  have  demonstrated  that  the  method  is  efficient  and  have 
confirmed  the  theory  that  it  is  also  stable. 

Our  new  approach  is  based  on  considering  a  form  of  the  3-D  wave  equation  as  an 
ordinary  differential  equation  with  respect  to  range.  Then  a  formal  expression  of  the 
solution  is  written  in  terms  of  the  exponential  of  the  square-root  of  some  operator.  The 
artifice  used  in  this  paper  is  to  approximate  this  exponential  in  a  clever  way  by  the  product 
of  two  rational  functions  of  the  type  (1,1).  As  a  consequence  the  resulting  ODE-integration 
process,  requires  only  two  successive  tridiagonal  system  solutions.  The  theory  shows  that 
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our  method  is  unconditionally  stable.  Moreover,  it  is  so  accurate  that  larger  step-sizes  can 
be  afforded  resulting  in  substantial  savings  in  computational  times.  This  has  been  widely 
confirmed  by  the  numerical  tests.  Moreover,  angles  of  propagation  as  wide  as  31  dgrees 
have  been  accurately  handled. 
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